/*
    Copyright 2006-2012 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Integrates the data cubes into broadband filter images. \ingroup broadband

/** \defgroup broadband Broadband color integration (broadband)

    The broadband executable takes the output from mcrx and integrates
    the data cubes into broadband color images.  It also integrates
    the integrated SEDs into integrated magnitudes.  

    The configuration options and output file format for broadband are
    described on the \ref broadband-readme page.
*/

#include "config.h"
#include "mcrx-types.h"
#include "blitz-fits.h"
#include "fits-utilities.h"
#include "interpolatort.h"
#include "preferences.h"
#include "hpm.h"
#include "counter.h"
#include "misc.h"
#include "constants.h"
#include "mcrx-debug.h"
#include "boost/lexical_cast.hpp"
#include "vecops.h"
#include "mcrxversion.h"

using namespace mcrx;
using namespace blitz;
using namespace CCfits;
using namespace boost;

typedef interpolator< T_float , T_float , 1> T_interpolator;


/** Calculates the effective width of the filters. \ingroup broadband */
array_1 calculate_ewidth_lambda (const array_2& passbands,
				 const array_1& lambda)
{
  // this is the effective width of the filter, 
  // Integral dlambda*R(lambda)
  return array_1 
    (integrate_quantities_log(passbands, lambda));
}

/** Calculates the effective width of the filters in nu
    units. \ingroup broadband */
array_1 calculate_ewidth_nu (const array_2& passbands,
			     const array_1& lambda)
{
  // c*Integral dlambda*R(lambda)/lambda^2
  return array_1
    (constants::c0*
     integrate_quantities_log(passbands (tensor::i, tensor::j)/
			      (lambda(tensor::j)*lambda(tensor::j)),
			      lambda));
}

/** Calculates the effective wavelength of the filters as Integral
    dlambda*lambda*R(lambda)/ewidth. \ingroup broadband */
array_1 calculate_elambda (const array_2& passbands,
			   const array_1& lambda)
{
  // this is the effective wavelength of the filter, 
  return array_1
    (integrate_quantities_log(passbands (tensor::i, tensor::j)*
			      lambda(tensor::j),
			      lambda)/
     calculate_ewidth_lambda(passbands, lambda));
}

/// Loads the list of filters.  \ingroup broadband
vector<string> load_filter_list (const Preferences& p)
{
  // Read filter list
  const string list_file_name = p.getValue("filter_list", string ());
  cout << "Reading filter bands from file" << list_file_name << endl;
  ifstream list_file (list_file_name .c_str()); 
  if (!list_file) {
    cerr << "Error: Could not open filter file" << endl;
    exit (1);
  }
#ifdef BROKEN_ISTREAM_ITERATOR
  istream_iterator<string> a(list_file);
  istream_iterator<string> b;
  return vector<string> (a,b);
#else
  return vector<string> (istream_iterator<string> (list_file),
		         (istream_iterator<string> ()));
#endif

}


/** Loads the filter files and interpolates the transmission for the
    wavelength points in the SED. \ingroup broadband */
array_2 load_filters (vector<string>& filter_names, const array_1& lambda,
		      Preferences& p)
{
  const T_float lambda_conv = 
    p.getValue("filter_lambda_conversion", T_float ());
  const string filter_directory = 
    word_expand(p.getValue("filter_file_directory", string ()))[0];
  
  array_2 filter_curves (filter_names.size(), lambda.size());

  // Read filter curves
  cout << "Reading transmission curves:" << endl;
  for (int i = 0; i <filter_names.size(); ++i) {
    const string filter_file_name =filter_directory + filter_names [i];
    cout << "  " << filter_file_name << endl;
    ifstream filter_file (filter_file_name.c_str());
    if (!filter_file) {
      cerr << "Error: Could not open filter file " << filter_file_name << endl;
      exit (1);
    }
    vector<T_float> wavelengths, transmission;
    T_float t1, t2;
    while (filter_file.peek() == '#') {
      string junk;
      getline (filter_file, junk);
    }
    while (filter_file >> t1 >> t2) {
      wavelengths.push_back(t1 );
      transmission.push_back(t2 );
    }
    assert (wavelengths.size() ==transmission.size());
    if (wavelengths.size() < 2) {
      cerr << "Error: Too few points in transmission curve" << endl;
      exit (1);
    }
    transform (wavelengths.begin(),wavelengths.end(),wavelengths.begin(), 
	       bind2nd(multiplies<T_float> (), lambda_conv));

    T_interpolator ti;
    ti.setAxis(0, wavelengths) ;
    for (int j = 0; j < transmission.size(); ++j)
      ti.setPoint(j, transmission [j]);

    for (int j = 0; j < lambda.size(); ++j)
      if ((lambda (j)  < ti.lower_range () [0]) ||
	  (lambda (j)  > ti.upper_range() [0]))
	filter_curves (i, j) = 0;
      else
	filter_curves (i, j) = ti.interpolate(lambda (j));
  }
  cout << "\nRead in filters\n";

  DEBUG(1,ofstream of("filters.txt");of << filter_curves.extent(firstDim) << '\t' << filter_curves.extent(secondDim) << endl; for(array_2::iterator i=filter_curves.begin(); i!=filter_curves.end();++i) of << *i << '\t'; of<< endl;);

  return filter_curves;
}


/** This function takes redshift and wavelength and returns Lyman
    series forest attenuation from Madau 1995. */
array_1 LyAlphaForest_attenuation(const array_1& lambda, double z)
{
  const int nlines = 4;
  // Wavelengths of the first four lyman lines
  const double Ly_lambda[nlines] = {121.6e-9, 102.6e-9, 97.3e-9, 95.0e-9}; 
  const double Ly_A[nlines]      = {0.0036, 1.7e-3, 1.2e-3, 9.3e-4};

  array_1 tau(lambda.shape());
  tau=0;
  const array_1 lambda_rest (lambda/(1+z));

  // Add up the optical depth for the lines
  // Looping through all the lines
  for(int i=0;i<nlines;++i)
    // If the wavelength in question is shorter than the wavelength of the line 
    // Then the total optical depth is increased by  . . . some amount
    tau += where(lambda_rest<Ly_lambda[i],
		 Ly_A[i]*pow(lambda/Ly_lambda[i],3.46),
		 0);

  // And now the Lyman continuum absorption

  // No wavelength is emitted at longer wavelength than this
  const double Ly_limit = 91.2e-9;  
  const array_1 xc (where (lambda/Ly_limit>1, lambda/Ly_limit, 1));
  const double xem = 1+z;

  tau+= where(xem>xc,
	      0.25*pow(xc,3.0)*(pow(xem,0.46)-pow(xc,0.46)) +
	      9.40*pow(xc,1.5)*(pow(xem,0.18)-pow(xc,0.18)) +
	      -0.70*pow(xc,3.0)*(pow(xc,-1.32)-pow(xem,-1.32)) +
	      -0.0023*(pow(xem,1.68)-pow(xc,1.68)),
	      0);

  return array_1(exp(-tau));
}



/** Integrates a data cube into images in the specified
    filters. \ingroup broadband */
void integrate_images (FITS& input, FITS& output, 
		       vector<string>& input_hdus, const string& output_HDU,
		       const array_2& passbands, const array_1& ewidth_lambda, 
		       const array_1& ewidth_nu, 
		       const array_1& delta_lambda,
		       const array_1& attenuation)
{
  Array<float, 3> image;

  for (vector<string>::iterator i = input_hdus.begin();
       i != input_hdus.end(); ++i) {
    cout << "Reading " << *i << endl;
    try {
      ExtHDU& camera_HDU = open_HDU (input, *i);
      Array<float, 3> temp_image;
      if (input_hdus.size() > 1) {
	read (camera_HDU, temp_image);
	if (image.size() == 0) {
	  image.resize(temp_image.shape() );
	  image = 0;
	}
	image += temp_image;
      }
    else
      read (camera_HDU, image);
    }
    catch (CCfits::FITS::NoSuchHDU&) {
      cerr << "Could not open HDU " << *i << ", skipping" << endl;
    }
  }
  
  if(image.size()==0) {
    cout << "No data found for output to " <<  output_HDU << endl;
    return;
  }

  cout << "Calculating passband fluxes"<<endl;

  firstIndex i;
  secondIndex j;
  thirdIndex k;
  fourthIndex l;

  // Calculate the f_lambda_effective for the filters.

  // f_lambda*delta_lambda === f_nu*delta_nu, so it has no effect when
  // doing this integration. The net result will always be that it
  // goes as (1+z)^-4, plus the effect of the bandpass width
  // shifting. It's the units of the equivalent width that determines
  // whether the result is in lambda or nu units, but that is an
  // observer quantity and cares not about redshift
  Array<float, 3> broadband_flux 
    (sum (image(i,j,l)*
	  attenuation(l)*
	  passbands(k,l)*
	  delta_lambda(l), l)/ewidth_lambda (k)); 

  // create CAMERAi-BROADBAND HDU
  ExtHDU* broadband_HDU;
  cout << "Writing " << output_HDU << endl;
  try {
    broadband_HDU = &open_HDU (output, output_HDU) ;
  }
  catch (CCfits::FITS:: NoSuchHDU&) {
    std::vector<long> naxes;
    naxes.push_back(broadband_flux.extent(firstDim ) );
    naxes.push_back(broadband_flux.extent(secondDim ) );
    naxes.push_back(broadband_flux.extent(thirdDim ) );
    broadband_HDU = output.addImage (output_HDU, FLOAT_IMG, naxes);
	
    broadband_HDU ->writeComment("This HDU contains the camera images integrated over broadband filters.  The third dimension is filter, see HDU BROADBAND for which filters.");
  }
  write (*broadband_HDU, broadband_flux);

  // To aid the user, we write the conversion factor from _lambda to
  // _nu as keywords, though this is more easily accessible by using
  // ewidth in the FILTERS HDU.
  for (int i=0; i<ewidth_lambda.size(); ++i) {
    // prevent NaN here if a filter is completely outside of the
    // simulated wavelengths so that both ewidths are 0.
    const T_float conv = (ewidth_lambda(i)==ewidth_nu(i))? 1.0 : 
      ewidth_lambda(i)/ewidth_nu(i);
    broadband_HDU->addKey("to_nu_factor" + lexical_cast<string>(i+1),
			  conv,
			  "Conversion factor to nu units for band " +
			  lexical_cast<string>(i+1));
  }

  // copy keywords (e.g. WCS, units, normalization) from the first input HDU
  ExtHDU& input_HDU = open_HDU(input, input_hdus[0]);
  input_HDU.readAllKeys();
  broadband_HDU->copyAllKeys(&input_HDU);
}


/** Integrates an SED into integrated magnitudes in the specified
    filters. \ingroup broadband */
void integrated_magnitudes (FITS& input, FITS& output, 
			    const string& input_col, 
			    const string& L_lambda_eff_col,
			    const string& abmag_col,
			    const array_2& passbands,
			    const array_1& ewidth_lambda,
			    const array_1& ewidth_nu,
			    const array_1& lambda,
			    const array_1& attenuation,
			    T_float abmag_zp)
{
  // We do this by using the SED's in the INTEGRATED_QUANTITIES HDU
  array_1 L_lambda;
  ExtHDU& iq_HDU = open_HDU (input, "INTEGRATED_QUANTITIES");
  try {
    read (iq_HDU.column(input_col), L_lambda);
  }
  catch (CCfits::Table::NoSuchColumn&) {
    cerr << "No spectrum for " << input_col << ", skipping" << endl;
    return;
  }


  // Calculate the effective luminosity integrated over the filter
  // bands. See the image integration for the units of this
  // expression. Specifically, we use the ewidth_lambda to get L_lambda_eff.
  const array_1 broadband_L_lambda_eff
    (integrate_quantities_log(
			      L_lambda(tensor::j)*	  
			      attenuation(tensor::j)* 
			      passbands(tensor::i, tensor::j),
			      lambda)/ewidth_lambda(tensor::i));
     
  ExtHDU& filters_HDU = open_HDU (output, "FILTERS");
  try { filters_HDU.column(L_lambda_eff_col );}
  catch (CCfits::Table::NoSuchColumn&) { 
    filters_HDU.addColumn(Tdouble, L_lambda_eff_col,
			  1, iq_HDU.column(input_col ).unit() );
  }
  write (filters_HDU.column(L_lambda_eff_col),
	 broadband_L_lambda_eff, 1);

  // and now calculate the absolute AB magnitude
  // We do this by calculating the flux assuming the object is at 10 pc.
  // We need to convert pc to Sunrise units (m); this is where the 85.19
  // comes from, as 85.19 = 2.5log10( 4pi (pc/m)^2 )
  const array_1 ab_mag
    (-2.5*log10 (broadband_L_lambda_eff*ewidth_lambda/(ewidth_nu*abmag_zp)) + 
     5 + 85.19);

  // and write
  try { filters_HDU.column(abmag_col );}
  catch (CCfits::Table::NoSuchColumn&) { 
    filters_HDU.addColumn(Tdouble, abmag_col, 1, "" );
  }
  write (filters_HDU.column(abmag_col), ab_mag, 1);
}


/** Writes a history blurb to the specified HDU. \ingroup broadband */
void write_history (HDU & output) 
{
  ostringstream history;
  history << "This file was created by Sunrise (broadband) version "
	  << mcrx::version() << ", on ";
  time_t now = time (0) ;
  char*now_string = ctime (& now);
  now_string [24] = 0;
  history << now_string << " by user ";
  const char* u=getenv ("USER");
  assert(u != 0);
  string user (u);
  const char* h=getenv ("HOST");
  if(h == 0)
    h = getenv ("HOSTNAME");
  assert(h != 0);
  string host (h);
  history << user << " on host "<< host << ".";
  output.writeHistory(history.str());
}


/** The main function for the broadband executable. \ingroup broadband */
int main (int argc, char**argv)
{
  cout << "broadband (Sunrise) version " << mcrx::version() << ".\nCopyright 2004-2012 Patrik Jonsson.\nThis program comes with ABSOLUTELY NO WARRANTY. Sunrise is free software,\nand you are welcome to redistribute it under certain conditions. See the\nGNU General Public License for details.\n\n";

  if (argc != 2) {
    cout << "Usage: broadband <config file>\n";
    exit (1);
  }
  cout << "\nBroadband Filter Integration Program\n";

  Preferences p;

  const string parameter_file (argv [1]);
  if (!p.readfile (parameter_file)) {
    cerr << "Error opening parameter file: " << parameter_file << endl;
    exit (1);
  }
  else
    cout << "Reading parameter file: " << parameter_file << endl;

  // check if hpm should be enabled (default is disabled)
  if (p.defined("use_hpm") && p.getValue("use_hpm", bool ()))
    hpm::enable ();
  else
    hpm::disable ();
  //hpm::hpminit(0, "broadband");

  // check if counters should be disabled (default is enabled)
  if (p.defined("use_counters") && !p.getValue("use_counters", bool ()))
    counter::enable = false;
  else
    counter::enable = true;

  // check if CCfits should be verbose (default is not)
  if (p.defined("CCfits_verbose"))
    CCfits::FITS::setVerboseMode (p.getValue("CCfits_verbose" , bool ()));
    
  // Get redshift from parameter file
  const T_float z = p.defined("redshift") ? 
    p.getValue("redshift", T_float()) : 0;
  cout << "Using redshift " << z << endl;

  p.setValue("pwd", string (getenv ("PWD")), "Working directory") ;
  p.setValue("configfile", parameter_file , "Configuration file") ;

  // check if we should calculate images or just integrated magnitudes
  const bool skip_images =
    p.defined("calculate_images") && !p.getValue("calculate_images", bool ());

  const string input_file (p.getValue("input_file", string ()));
  const string output_file (strip_suffix (p.getValue("output_file", string ()))
			    + ".fits");
  // Load filter list
  vector<string> filter_names = load_filter_list (p);
  if (filter_names.empty()) {
    cerr << "No filter bands defined, nothing to do!" << endl;
    exit (0);
  }

  auto_ptr<FITS> input (new FITS (string (input_file), Read));    

  // Read wavelengths from input file
  cout << "Reading wavelengths" << endl;
  ExtHDU& lambda_hdu = open_HDU (*input, "LAMBDA");
  const string lambda_unit = lambda_hdu.column("lambda" ).unit();
  array_1 lambda;
  read (lambda_hdu.column("lambda" ), lambda);
  lambda = lambda*(1+z);

  // *** NOTE THAT LAMBDA NOW IS *OBSERVED* LAMBDA, NOT REST ***

  const array_1 delta_lambda = delta_quantity (lambda);
  
  // load the filter passbands and return an array of (filter,
  // wavelength)
  const array_2 passbands = load_filters (filter_names, lambda, p);
  // calculate effective width and write to file
  const array_1 ewidth_nu= calculate_ewidth_nu (passbands, lambda);
  const array_1 ewidth_lambda= calculate_ewidth_lambda (passbands, 
							lambda);

  // Deal with units
  ExtHDU& iq_HDU = open_HDU (*input, "INTEGRATED_QUANTITIES");
  const string L_lambda_unit = iq_HDU.column("L_lambda").unit(); 
  string L_nu_unit = L_lambda_unit;
  const string nu_unit = "Hz";
  string::size_type pos = L_lambda_unit.find ("/" + lambda_unit);
  if (pos == string::npos) {
    // OK, now this is weird.  These units SHOULD be consistent.
    cerr << "Warning: wavelength and luminosity units are not consistent: " << lambda_unit << " " << L_lambda_unit << endl;
    L_nu_unit += lambda_unit;
  }
  else
    L_nu_unit = L_nu_unit.replace(pos, 1+ lambda_unit.size(),
				  string ("") );
  string conv_unit;
  if (lambda_unit == "m") {
    L_nu_unit += "/Hz";
    conv_unit = "m/Hz";
  }
  else {
    // OK, problem.  We don't know how to deal with this, so just make it consistent
    L_nu_unit += "/Hz/m*" + lambda_unit;
    conv_unit = lambda_unit + "*" + lambda_unit + "/m/Hz";
  } 
  T_float abmag_zp = 1;
  if (L_nu_unit == "W/Hz")
    abmag_zp = 3.63e-23; // The AB zero point is 3.63e-23 W/m2Hz
  else
    cerr << "Warning: Don't know AB magnitude zero point in units "
	 << L_nu_unit << ".  Using 1, probably won't make sense!"  <<  endl;

  // Calculate an attenuation factor which is used to attenuate the
  // flux from the objects. This can include lyman alpha forest
  // absorption as well as surface brightness dimming.
  array_1 attenuation (lambda.shape());
  array_1 attenuation_integrated (lambda.shape());

  // The correct power of 1+z here is debated. I now claim that for I_lambda it is 5!
  // * time dilation of photon arrival
  // * stretching of dlambda spreading out photons in wavelength
  // * redshifting of photon energy
  // * 2 powers from angular diameter distance
  // Note that one of those powers will be negated if we integrate
  // over the entire wavelength range, ending up with the canonical 4
  // for bolometric surface brightness.
  attenuation = 1.0*pow(1+z, -5.);

  // For the integrated magnitudes, the effect is only one power of
  // 1+z which must be applied to counteract the fact that we are
  // redshifting the spectrum while keeping the total energy constant.
  attenuation_integrated = 1.0*pow(1+z, -1.);

  // If we are using lyman absorption, calculate that now
  // LyA forest attenuation only depends on wavelengths and z so we
  // can calculate the attenuation factor now once and for all
  if (p.defined("use_lyman_absorption") ? 
      p.getValue("use_lyman_absorption", bool()) : false) {
    cout << "Including Lyman Alpha forest absorption in calculation." << endl;
    attenuation *= LyAlphaForest_attenuation(lambda, z);
    /// \todo no forest for integrated magnitudes?
    //attenuation_integrated *= LyAlphaForest_attenuation(lambda_integrated, z);
  }

  // now open the output file
  // see if we are supposed to resume 
  auto_ptr<FITS> output;

  // check if previously started and should be resumed, or if done
  int resume_camera = 0;
  // first try to open an existing file
  try {
    output.reset(new FITS (output_file, Read));

    // if that succeeds, see if we are done
    bool done=false;
    try {
      output->pHDU ().readKey("BROADBAND_COMPLETE", done);
    }
    catch (...) {}
    if (done) {
      cout << "\nBroadband already complete" << endl;
      return 0;
    }
    
    // if we are not done, see if we should resume
    try {
      output->pHDU ().readKey("BBRESUME", resume_camera);
    }
    catch (...) {}
    if (resume_camera > 0) {
      cout << "\nResuming with camera " << resume_camera << endl;
    }
  }
  catch (...) {
    // failed opening input file, that means we should create it
    cout << "Creating output file " << output_file << endl;
    // "copy primary HDU" constructor
    FITS output ("!" + output_file, *input);
    // Write the primary keywords
    output.pHDU().addKey("FILETYPE", string ("BROADBAND"),
			 "File contains broadband colors");
    output.pHDU().addKey("DATATYPE", string ("GRID"),
			 "Data is on a grid");
    output.pHDU().addKey("FILTERCODE", string ("BROADBAND"),
			 "See HDU BROADBAND for filter setup");

    // output keywords to HDU "broadband"
    ExtHDU* broadband_HDU = output.addTable("BROADBAND", 0);
    write_history (*broadband_HDU);
    broadband_HDU->writeDate();
    p.write_fits(*broadband_HDU );
    
    // copy the other HDU's from the input file
    cout << "Copying HDU's from input file" << endl;
    try {  
      const CCfits::ExtMap extensions = 
	const_cast<const FITS*>(input.get())->extension();  
      for (multimap< string, ExtHDU*>::const_iterator i = extensions.begin();
	   i != extensions.end(); ++i) {
	// exclude CAMERAi (-SCATTER/-NONSCATTER/-ATTENUATION/etc) HDU's
	const string ext_name (i->first);
	if (ext_name.find(string ("CAMERA")) != string::npos &&
	    ext_name.find(string ("-PARAMETERS")) == string::npos)
	  continue;
	// also exclude these:
	if ((ext_name == "GRIDSTRUCTURE") ||
	    (ext_name == "LAMBDA") ||
	    (ext_name == "DEPOSITION") ||
	    (ext_name == "INTENSITY") ||
	    (ext_name == "DUSTINTENSITY") ||
	    (ext_name == "CELLSEDS") ||
	    (ext_name == "RANDOM_STATE"))
	  continue;
	cout << "  " << ext_name << endl;
      
	output.copy(input->extension(ext_name) );
      }
    }
    catch(...) {
      cerr << "Error copying HDU's, output file is incomplete!"  << endl;
    }

    // Create FILTERS HDU and write filter bands
    ExtHDU* filters_HDU = output.addTable("FILTERS", 0,
					  vector<string> (), vector<string> (),
					  vector<string> ());
    filters_HDU->addColumn(Tstring, "filter", 30, "Filter name" );
    filters_HDU->column("filter" ).write(filter_names, 1);

    filters_HDU->addColumn(Tdouble, "lambda_eff", 1, lambda_unit );
    write (filters_HDU->column("lambda_eff"), 
	   calculate_elambda(passbands, lambda), 1);
    filters_HDU->addColumn(Tdouble, "ewidth_lambda", 1, lambda_unit );
    write (filters_HDU->column("ewidth_lambda"), 
	   calculate_ewidth_lambda(passbands, lambda), 1);
    filters_HDU->addColumn(Tdouble, "ewidth_nu", 1, nu_unit );
    write (filters_HDU->column("ewidth_nu"), ewidth_nu, 1);
  } // initialize output file

  int n_cameras;
  output.reset(0 );
  output.reset(new FITS (output_file, Write));
  open_HDU (*output, "MCRX").readKey("N_CAMERA", n_cameras);

  cout << "\nIntegrating scatter images" << endl;
  // *** LOOP OVER SCATTER CAMERAS ***
  if (resume_camera < n_cameras) {
    for (int i = resume_camera; i < n_cameras; ++i) {
      std::ostringstream ost;
      ost << i;
      const string camera_numstr=ost.str();

      output.reset(0 );
      output.reset(new FITS (output_file, Write));

      // integrate dust images
      input.reset (0);
      input.reset (new FITS (string (input_file), Read));

      vector<string> hdu_names;
      hdu_names.push_back("CAMERA" + camera_numstr);
      hdu_names.push_back("CAMERA" + camera_numstr + "-IR");
      if (!skip_images)
	integrate_images (*input, *output, hdu_names,
			  "CAMERA" + camera_numstr+ "-BROADBAND",
			  passbands,
			  ewidth_lambda, ewidth_nu,
			  delta_lambda,
			  attenuation);
		      
      // Calculate integrated passband flux as well
      cout << "Calculating integrated magnitudes" << endl;

      integrated_magnitudes(*input, *output, 
			    "L_lambda_out"+camera_numstr,
			    "L_lambda_eff" + camera_numstr,
			    "AB_mag" + camera_numstr,
			    passbands,
			    ewidth_lambda, 
			    ewidth_nu, 
			    lambda,
			    attenuation_integrated,
			    abmag_zp);

      // in case program is killed, write where we were to the keyword
      // so we can resume
      output->pHDU().addKey("BBRESUME", i+ 1,
			    "Interrupted at this camera" );

    }
    resume_camera = 0;
  }
  else
    // if we are resuming nonscatter images, adjust the value of the
    // camera
    resume_camera -= n_cameras;
  
  cout << "Integrating nonscatter images" << endl;
  // *** LOOP OVER NONSCATTER CAMERAS ***
  for (int i = resume_camera; i < n_cameras; ++i) {
    std::ostringstream ost;
    ost << i;
    const string camera_numstr=ost.str();

    output.reset(0 );
    output.reset(new FITS (output_file, Write));
    
    // and then no scatter images
    input.reset (0);
    input.reset (new FITS (string (input_file), Read));

    vector<string> hdu_names;
    hdu_names.push_back("CAMERA" + camera_numstr + "-NONSCATTER");
    if (!skip_images)
      integrate_images (*input, *output,
			hdu_names,
			"CAMERA" + camera_numstr+ "-BROADBAND-NONSCATTER",
			passbands, ewidth_lambda, ewidth_nu,
			delta_lambda, 
			attenuation);
    
    // Calculate integrated passband flux as well
    cout << "Calculating integrated magnitudes" << endl;

    // The nonscatter SED has the same wavelengths as the images... (hoaky)
    integrated_magnitudes(*input, *output, 
			  "L_lambda_nonscatter"+camera_numstr,
			  "L_lambda_eff_nonscatter" + camera_numstr,
			  "AB_mag_nonscatter" + camera_numstr,
			  passbands,
			  ewidth_lambda, 
			  ewidth_nu, 
			  lambda,
			  attenuation_integrated,
			  abmag_zp);

    // in case program is killed, write where we were to the keyword
    // so we can resume
    output->pHDU().addKey("BBRESUME", n_cameras + i+ 1,
			  "Interrupted at this camera" );
  }

  output->pHDU().deleteKey("BBRESUME");
  output->pHDU().addKey("BROADBAND_COMPLETE", true,
			"Broadband calculation completed.");

}
